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ABSTRACT 

We present a novel numerical implementation of radiative transfer in the cosmological 
smoothed particle hydrodynamics (SPH) simulation code GADGET. It is based on a 
fast, robust and photon-conserving integration scheme where the radiation transport 
problem is approximated in terms of moments of the transfer equation and by using a 
variable Eddington tensor as a closure relation, following the 'OTVET'-suggestion of 
Gnedin & Abel. We derive a suitable anisotropic diffusion operator for use in the SPH 
discretization of the local photon transport, and we combine this with an implicit solver 
that guarantees robustness and photon conservation. This entails a matrix inversion 
problem of a huge, sparsely populated matrix that is distributed in memory in our 
parallel code. We solve this task iteratively with a conjugate gradient scheme. Finally, 
to model photon sink processes we consider ionisation and recombination processes 
of hydrogen, which is represented with a chemical network that is evolved with an 
implicit time integration scheme. We present several tests of our implementation, 
including single and multiple sources in static uniform density fields with and without 
temperature evolution, shadowing by a dense clump, and multiple sources in a static 
cosmological density field. All tests agree quite well with analytical computations or 
with predictions from other radiative transfer codes, except for shadowing. However, 
unlike most other radiative transfer codes presently in use for studying reionisation, 
our new method can be used on-the-fly during dynamical cosmological simulation, 
allowing simultaneous treatments of galaxy formation and the reionisation process of 
the Universe. 
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1 INTRODUCTION 



The absence of Gunn- Peterson troughs in the spectra of high 
redshift quasars up to ^ ^ 6 (jPan et al.l l2006l : I White et al.l 
l2003l ) suggests that hydrogen is highly ionised at low red- 
shift, with a volume av eraged neutral hydrogen fraction of 
xm ^ 10~^ (|Fanll2007l V The current generation of simula- 
tion codes for cosmological structure formation calculates 
the self- gravity of dark matter and cosmic gas, and the fluid 
dynamics of the cosmic gas, but radiation processes are typi- 
cally not taken into account, or only at the level of a spatially 
uniform, externally imposed background field. However, we 
know that the radiation field has been highly inhomoge- 
neous during certain phases of the growth of structure, and 
may have in fact provided i mportant feedback effects for 
galaxy form ation (e.g. illiev et al.ll2005l : lYoshida et akl lioOTl : 
ICroft Altav 2008 ). Therefore, it would be ideal to be able 
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to calculate self- consistent simulations that simultaneously 
follow cosmic structure growth and the ionisation history of 
the Universe. 

However, the high calculation cost and complexity of ra- 
diative transfer has so far prevent ed such simulations, apart 
from very few exceptions (e.g. iGnedin Ostrikerl Il997l : 
iKohler et al.ll2007l : IShin et al.l l2008l : IWise k Abelll2008l ) . In- 
stead, the cosmic reionisation process has been most often 
treated separately in pos tprocessing, based on static simu- 
lated density fields (e.g 
20041 : Illiev et al.1 I2OO6I ; 



Ciardi et ai] |2003l : ISokasian et al 



Zahn et al.1 120071 : ICroft Alt~ 



Li et al.l I2OO8I I . It is the goal of this study to de- 
velop a new numerical scheme that makes self- consistent 
radiation-hydrodynamic simulations possible based on the 
high- resolution Lagrangian cosmol ogical code GADGET-3 
(jSpringel et al.ll200ll : ISDringelll2005l l. 

A large number of different approximate methods have 
been developed to make problems of radiative transfer 
numerically tractable, which are in their full generality only 
solvable analytically for very few special cases. The known 
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numerical methods include lon^ chara ct eristi cs ray-tracing 
schemes (iMihalas &: Weibel MihalasI Il984l: lAbel et al 



1998 



20021 



f 



ISokasian et al. 2001; Cen 
Kazoumov & Cardall 2005i) , 



ray-t ra cing schemes (iKunasz Aueu 



20011: iMellema et al 
Alvarez et al l l2006l: 



120021 : lAbel Wandelt 
shor t characteristic 
Nakamoto et alJ 



IWhalen &: Norman 



Altav et al 



Ciardi et al.1 l200ll: iMaselli et al l 12003). hybrid ray-tracing 
schemes (jmikhorst et al.1 l2006l : LTrac Cen 2007) , mo- 



ment methods and direct solvers (pMellema et al l Il998l: 
Gnedin Abelll200ll: IShaDiro et al."20 04l: IWhitehouse et al.1 
20051 : lAubert k Tevss ier 2008; Finlat or et al.1 l2008l ) and 



other particle-b ased tr ansport methods (|Ritzerveld et all 
l2003l : [Pawlik k Schavd lioOS). 

In the long characteristics method each test cell in the 
computational volume is connected to all other relevant 
cells. Then the equation of radiative transfer is integrated 
individually from that test cell to each of the selected cells. 
While this method is relatively simple and straight-forward, 
it is also very time consuming, since it requires N'^ inter- 
actions between the cells. Moreover, parallelization of this 
approach is cumbersome and requires large amounts of data 
exchange between the different processors. Short character- 
istic methods integrate the equation of radiative transfer 
along lines that connect nearby cells. Here a cell is connected 
only to its neighboring cells, and not to all other cells in the 
computational domain. This reduces the redundancy of the 
computations and makes the scheme easier to parallelize. 

Stochastic integration methods, specifically Monte 
Carlo methods, employ a ray-tracing strategy where the rays 
are discretized into photon packets. For each photon packet, 
the emission, its frequency and its direction of propagation 
are determined by sampling the appropriate distribution 
functions that have been assigned in the initial conditions. A 
particular advantage of this approach is that comparatively 
few approximations to the radiative transfer equations need 
to be made, so that the quality of the results is primarily a 
function of the number of photon packets employed, which 
can be made larger in proportion to the CPU time spent. 
A disadvantage of these scheme is the comparatively high 
computational cost and the sizable level of noise in the dis- 
cretized radiation field. 

Using moments of the radiative transfer equations in- 
stead of the full radiative transfer equation can lead to very 
substantial simplifications that can drastically speed up the 
calculations. In this approach, the radiation is represented 
by its mean intensity field throughout the computational do- 
main. Instead of following rays, the moment equations are 
solved directly on the grid or for the particles. Due to its 
local nature, the moment approach is comparatively easy to 
parallelize, but it of course depends on the nature of the 
investigated problem whether the simplifications made still 
provide sufficient accuracy. 

In this paper we develop a moments method that 
is closely related in spirit to the Optically Thin Vari- 
able Eddington Ten sor (OTVET) scheme proposed by 
iGnedin Abel (|2001 ). However, we try to implement it di- 
rectly on top of the irregular set of positions sampled by the 
particles of Smoothed Particle Hydrodynamics (SPH) sim- 
ulations, and we use quite different numerical techniques to 
solve the resulting transport equations. In OTVET, the sys- 
tem of moment equations is closed by estimating the local 



Eddington tensor with a simple optical thin approximation, 
i.e. one pretends that all sources of light are 'visible' at a 
given location. Once the Eddington tensors are found, the 
local radiation transfer reduces to an anisotropic diffusion 
problem. The particular attraction of this moment-based 
formulation is that it is potentially very fast, allowing a di- 
rect coupling with cosmological hydrodynamic simulations. 
In particular, if a rapid method for calculating the Edding- 
ton tensors can be found, the scheme should be able to easily 
deal with an arbitrary number of sources. Also, the radia- 
tion intensity field does not suffer from the Poisson shot 
noise inherent in Monte Carlo approaches. Together with 
the local nature of the diffusion problem, this makes this 
approach particularly attractive for trying to address the 
cosmological reionisation problem with self-consistent simu- 
lations of galaxy formation, since it is likely that low-mass 
star-forming galaxies of high number density play an impor- 
tant role for the reionisatio n process. We therefor e adopt in 
this paper the suggestion o f iGnedin Abel (|200lh and work 
out an implementation of the OTVET scheme in SPH. As 
we shall see, this entails a number of numerical challenges in 
practice. We will describe our solutions for these problems, 
and carry out a number of tests to evaluate the accuracy of 
the resulting implementation. 

The near complete lack of analytical results for non- 
trivial radiative transfer problems makes it actually hard 
to vaUdate different numerical techniques and to compare 
their performance with each other. A very useful help in this 
respect is provided by the cosmolog!;ical radiative t ransfe r 
code comparison project, carried out by llliev et al.1 (|2006l ). 
In their paper, they present a comparison of 11 independent 
cosmological radiative transfer codes when applied to a va- 
riety of different test problems. A number of our tests are 
based on this study, which hence allows a comparison with 
results of these other codes. 

We start this paper with a brief introduction to the ra- 
diative transfer equations in Section [2] We then describe in 
Section [3] the moment-based method that is the basis for our 
approximate treatment of the radiation transfer problem. In 
Section m we elaborate in detail the numerical implementa- 
tion of this scheme in SPH. This is followed by a presentation 
of results for various test problems in Section [5] Finally, we 
conclude with a summary and an outlook in Section [6] In an 
appendix, we give for reference further details on our rate 
equations, and on the conjugate gradient approach. 



2 THE RADIATIVE TRANSFER PROBLEM 

2.1 Equation of radiative transfer 

Let us briefly derive the radiative transfer (RT) equation 
in comoving coordinates, which is also useful for introduc- 
ing our notation. Let /^(t,x,p) be the photon distribution 
function for comoving coordinates x and comoving photon 
momentum 

P — ' 

where a = a(t) is the cosmological scale factor, h is the 
Planck constant, u is the frequency of the photons, and fi is 
the unit vector in the direction of photon propagation. Then 
the number of photons in some part of the Universe is 
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dxdp/^(t,x,p) . 



(2) 



We can further define the phase- space continuity equa- 
tion for the distribution function = fj{t, x, p) of photons 

as 



d 



dt 



dt 



(3) 



Here the source and sink terms on the right hand side 
of the equation represent photon emission and absorption 
processes, respectively. We define the specific radiation in- 
tensity lu as the energy of photons in a frequency bin 
that pass through an area AA and sohd angle AQ for a time 
At. The specific intensity lu is then related to the photon 
distribution as follows 



(4) 



^^diyd^dAdt 
Substituting into equation JS]), rearranging and adding 
the proper absorption and emission terms, one obtains the 
following radiative transf er equation in comoving variables 
(jGnedin QstrikejlToOTl ) : 

c at a oyi c \ ov J 

where H = d/a is the Hubble expansion rate, is the 
absorption coefficient and ji, is the emission coefficient. 

Unfortunately, this full radiative transfer (RT) equation 
is in practice very difficult to solve in full generality. In par- 
ticular, the high dimensionality (comprised of 3 spatial vari- 
ables, 2 directional angles, 1 frequency variable, and time) 
of this partial differential equation makes a direct discretiza- 
tion on a mesh highly problematic. We will hence apply in 
Section [3] simplifications to the RT equation that yield an 
approximation that can be more easily calculated in cosmo- 
logical codes. 



2.2 Basic physics of hydrogen photoionisation 

If we consider pure hydrogen gas, the rate equations de- 
scribing photoionisation and photoheating processes become 
comparatively simple. For the most part we will restrict our- 
selves to this chemical composition in this paper, but we note 
that an extension of our formalism to include other elements 
is readily possible. In fact, we have already implemented he- 
lium as well, but we omit an explicit discussion of it in the 
following for the sake of simplicity. 

The photoionisation rate hon of hydrogen {H -\- hu ^ 

-\- e~) is given by 



k 



diy- 



hv 



(6) 



where hvo — 13.6 eV is the hydrogen ionisation potential 
and is the photoionisation cross-section: 

zy \4 exp{4 - [(4 tan"^e)/e]} 



■ exp(— 27r/e) 



for ^ z/o, 



(7) 



where (To — 6.30 x lO""*^^ cm^ and e — y^^^V^^o) — 1- 

The corresponding photoheating rate of hydrogen is 
given by 



r = riHi / dQ / diy^^jp^{hiy 



(8) 



where nni is the number density of neutral hydrogen. Fur- 
thermore, the change in the neutral gas density due to re- 
combinations is given by 

dnm 



dt 



aUeflp, 



(9) 



where a is the temperature-dependent recombination coeffi- 
cient. He is the electron number density and Up = nmi is the 
proton number density, which is in turn equal to the ionised 
hydrogen number density (for a pure hydrogen gas). 

The change of the density of the neutral gas due to 
ionisations is given by 

dnm .-.^x 
- -caoumn^, (10) 

where c is the speed of light and rij is the number density of 
ionising photons. Thus, the total change in the neutral gas 
density, due to ionisations or recombinations, is given by 

dnm 



dt 



aUeUmi — caoUmn^ 



(11) 



For all our calculations described in this 



pape 



use t he on-the-spot approximation ([Osterbrock &; Ferlandl 
I2OO6I ). i.e. photons emitted due to recombinations to ex- 
cited levels are re-absorbed immediately by neutral hydro- 
gen atoms in the vicinity. This behavior is described by the 
so called case-B recombination coefficient as- 



3 THE VARIABLE EDDINGTON TENSOR 
FORMALISM 

We now turn to a description of the moment-based approx- 
imation to the radiation transfer problem that we use in 
this study. The first three moments of the specific intensity, 
the mean intensity J^^ the radiation flux vector Fl^ and the 
radiation pressure tensor Pl\ are defined as follows: 



I = - I 

4n J 

^ - ^ [ 



dQI^, 



dQn^Iu, 



dQnn^Iiy 



(12) 



(13) 



(14) 



where n is a direction vector and the indices i and j run 
through the three elements of the vector in Cartesian space. 
We can further define /i*^, the so-called Eddington tensor, 
based on P^^ = Jjyh''^ . 

We can for the moment ignore the frequency derivative 
in the RT equation if we can assume that the Universe does 
not expand significantly before a photon is absorbed. With 
this simplification, the first moments of the RT equation 
take the form: 



ldJ^^ ^ IdFZ 
c dt a dx'^ 
1 dFi ^ IdJuh^ 
c dt a dx^ 
where 



(15) 
(16) 

(17) 



In the second moment equation ()16|) , we can ignore the term 
of the order and solve for the flux 
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= - — - 



1 IdJ^h'' 



(18) 



which we then insert back into equation (|15|) . This leads to 
the following approximation to the RT equation: 



dt 



c d 



]_djjiy_ 



(19) 



This form of the RT equation is already much simpler 
than the fully general form of equation (|5]). In particular, 
each of the terms in equation (|19|) has a simple physical 
interpretation. The time evolution of the local mean radi- 
ation intensity is given by a transport term, described by 
the anisotropic diffusion term on the right hand side, a sink 
term describing absorptions, and an emission term that ac- 
counts for sources. However, in order to be able to solve 
this equation an expression for the Eddington tensor K^^ is 
needed, which is left undefined by these moment equations. 
We therefore need to assume a certain form for the Edding- 
ton tensor, or in other words, a closure relation. 

F or the closure relation, we follow Gn edin k, Abell 
(|200lh and estimate the local Eddington tensor with an op- 
tically thin approximation. This means that we assume that 
a reasonable approximation to the Eddington tensor can be 
obtained by approximating all lines-of-sight to the sources 
as being optically thin. The radiation intensity pressure ten- 
sor P^^ in this optically thin regime can then be computed 
as 

-x')i(x-x')^- 



/ 



(x- 



and thus the Eddington tensor is given by 



TV(P)- 



(20) 



(21) 



Note that the Eddington tensor only determines in which 
direction the local radiation propagates, but the magnitude 
of the radiation intensity tensor is unimportant as far as 
the Eddington tensor is concerned. This means that even in 
situations where the lines-of-sight to the sources are not op- 
tically thin at all, one will often end up with fairly accurate 
estimates of the Eddington tensor based on equations (|20)) 
and (|2ip , simply because the radiation will typically mainly 
propagate away from the sources, even in optically thick 
cases. In particular, note that the above approximation is 
always correct for a single source. When there are multiple 
sources of equal strength, the optically thin approximation 
will weight the sources that are closest most strongly, in ac- 
cordance with the decay of the intensity. While this 
can be expected to result in reasonably accurate estimates 
of the Eddington tensor in many situations (especially in the 
vicinity of a dominating source), errors can certainly arise 
in particular situations, for example at locations that are 
equidistant from two sources of equal strength. How serious 
these errors are in problems of interest needs to be analyzed 
with appropriate test problems. 

As we describe later in more detail in Section 14. 8[ we 
note that equation (|2Qp can be accurately calculated with a 
hierarchical multipole approach similar to the one applied 
in gravitational tree algorithms. This allows a fairly efficient 
treatment of an arbitrarily large number of sources, which 
is a distinctive advantage of the moments based approach 
compared with other methods. 



3.1 Choice of convenient Lagrangian variables 

We will now rewrite the RT equations into a form that is 
more convenient for use with a Lagrangian method such as 
SPH. In particular, it is advantageous to pick variables in 
the numerical scheme that are normalized to unit mass, not 
unit volume. For example, if we express the ionisation state 
of the gas as the number density of ionised hydrogen per 
unit volume, then we have to readjust this number some- 
how any time we reestimate the local gas density (which 
may change if the gas moves around), otherwise the ionised 
fraction would change. However, if we use convenient vari- 
ables that are normalized to unit mass, we do not need to 
worry about such corrections. 

For chemical networks of hydrogen, it is convenient (and 
often done in practice) to express abundances relative to the 
total abundance of hydrogen nuclei: 



(22) 



Here — 0.76 is the cosmological mass fraction of hy- 
drogen and mp is the proton mass. In the following, we use 
the notation nni for neutral hydrogen, and nnii for ionised 
hydrogen, such that 



(23) 



In our actual numerical code, we will use a variable nnii 
to express the abundance of ionised hydrogen, defined as 



riHii ■ 



(24) 



where nnii is the ordinary number density of HII atoms (i.e. 
number of protons per unit volume) . Note that this quantity 
is now normalized to unit mass, as desired. In addition, it is 
dimensionless, which avoids numerical problems due to large 
numbers if we use astronomical length units. 

A similar reasoning also applies to the radiation inten- 
sity itself. In principle, the fundamental quantity we work 
with is the frequency dependent, angle- averaged mean in- 
tensity. However, we cannot afford to carry around a full 
spectrum with each fluid element in a hydrodynamical code. 
This would be too cumbersome and also is not really nec- 
essary of we are interested only in the reionisation problem. 
Instead, it is sufficient to store the intensity integrated over 
a narrow frequency interval around the ionisation potential 
of hydrogen. Or in other words, a more convenient quantity 
to work with would be something like the number density 
of photons capable of ionising hydrogen. We now formulate 
the relevant equations using this concept. 

In general, the photon number density is 



47rJ^ 



dzy. 



(25) 



However, we will only consider the spectrum in a small band 
around the frequencies of interest. For simplicity, we assume 
that the spectrum has the form 



(26) 



around the ionisation frequency z/q, where huo = 13.6 eV is 
the hydrogen ionisation potential. This form of the radiation 
intensity limits the spectrum to effectively just the hydrogen 
ionisation frequency. We therefore obtain this simple form 



Radiative transfer in GADGET 5 



1 47rJo 



(27) 



c huo 

for the number density of ionising photons. 

We also note that the absorption coefficient for ion- 
isation in the equation for radiation transport is 



(28) 



where cTi, is the cross-section for hydrogen ionisation. If we 
multiply the loss term HiuJu in the RT equation by An / (c hv) 
and integrate over i/, we get the so-called ionisation rate /cion, 
given by 



/ 47r J, 



(29) 



For our narrow spectrum, this leads to the simple expression 



A 4^ 7 

az/— — KvJi, — cToriHi , 



(30) 



where ctq is the cross section at the resonance. Another con- 
sequence of these definitions is that we can write the number 
density evolution of ionised hydrogen due to new reionisa- 
tions as 

driHii 



dt 



= c (To ^7 nui 



(31) 



The photon field loses energy at the same rate, i.e. the loss 
term for the radiation field should be of the form 



-cfjo riHi, 



(32) 



dn^ 
"dT 

which is also what the loss term in the RT equation gives in 
this notation. 

The above suggests that we can cast the moment-based 
RT equation into a more convenient form if we multiply it 
through with An / [chv) and integrate over v. This gives: 



first term on the right hand side is a diffusion like equa- 
tion, which is conservative, i.e. it leaves the total number 
of photons unchanged. The second term describes photon 
losses, and each photon lost will cause one hydrogen atom 
to be ionised. Finally, the third term is the source term, and 
describes the injection of new photons. This suggests that 
the total number of ionisations must always be equal to to 
the total number of photons lost. If we can maintain nu- 
merically accurate photon conservation, then this property 
should ensure a proper speed of the ionisation front even 
for relatively inaccurate time stepping, as the propagation 
of the front should largely be determined by the injection 
rate of photons at the source. 

The above suggests a simple possibility for treating 
the time evolution of the photon number density in each 
timestep in terms of three parts, corresponding to an 
operator- splitting^ or fractional- step approach: One may first 
inject new photons according to the source function, then 
transport photons conservatively by treating the diffusion 
part, and finally, advance the "chemical network" fean. [35)) 
by treating only ionisations and recombinations, making 
sure again that we do not lose any photons. The chemical 
equations can be easily "subcycled" or treated with an inte- 
grator for stiff differential equations, if needed, because they 
are completely local. On the other hand, the most expensive 
part of the time advance is given by the diffusion part. This 
not only involves a coupling with neighboring fluid elements 
but also cannot easily be integrated with an explicit time 
integration scheme, because the diffusion equation becomes 
easily unstable in this case. We will therefore treat this part 
with an implicit method. While involving an expensive it- 
eration scheme, this provides good stability and allows for 
comparatively large timesteps. 



1 dn^W' 
K dxi 



■ CKUj + S 



(33) 



where k, = aoum and the cosmological scale factor has been 
dropped for simplicity. The source function gives the rate 
per unit volume at which new ionising photons are produced. 

Finally, we also want to express the photon number 
density in terms of the local density of hydrogen atoms. 



n^/nti- Then we get: 



dxj 



K dxi 



■ C Ki Tl'-y ~\~ S' 



7- 



(34) 



This is the formulation of the RT equation that we imple- 
mented in this work in the simulation code GADGET-3. 

We have to augment equation (|34)) with the changes of 
the different chemical species as a result of interactions with 
the radiation field. If we consider only hydrogen, this is just: 



dt 



= caonu nuin 



7 ■ 



a nu rienmi 



(35) 



where a is the temperature-dependent recombination coef- 
ficient. If we have only hydrogen, we can set fie = nnii and 
rim = 1 — nmi- 

For hydrogen reionisation problems, we want to solve 
the two basic equations (|34|) and (|35p as efficiently, accu- 
rately and robustly as possible. We recall that the three 
terms of (|34|) have a straightforward interpretation. The 



4 NUMERICAL IMPLEMENTATION 

In this section we describe the numerical formalism we have 
implemented in order to solve the moment-based RT equa- 
tions coupled to the parallel Tree/SPH code GADGET-3, 
which is a significantly evolved and extended version of the 
pubUc GADGET-2 code (Springel 2005). We first give a very 
brief overview of the basic concepts of smoothed particle 
hydrodynamics (SPH), and then present a derivation of a 
new anisotropic diffusion operator in SPH, which is needed 
for the radiation transfer in moment form when a spatially 
varying Eddington tensor is used. We also explain how the 
diffusion equation can be integrated robustly in time based 
on an implicit scheme with an iterative sparse matrix solver. 
The time integration of the rate equations for ionisation and 
recombination also requires special methods because they 
involve stiff differential equations. Finally, we describe the 
calculation of the Eddington tensors, and how this can be 
best combined existing algorithms in the GADGET-3 code. 



4.1 The basics of smoothed particle 
hydrodynamics 

SPH is a widely used Lagrangian scheme that follows the 
evolution of g as properties bas ed on discrete tracer parti- 
cles (see, e.g., iMonaghanl 1 19921 , for a review). The particle 
properties are averaged, 'smoothed', over a kernel function. 
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yielding so-called kernel interpolants for the fluid properties 
based on a few sampling points. For example, the kernel- 
interpolant of a property {A) is given by 



j dr'A(r')H^(|r' -r|,/i), 



(36) 



where h is called the smoothing length and is defined such 
that the kernel W drops to zero for |r^| > /i. In a discretized 
form this equation becomes 



(^(rO) = V^A(r,)VK(|r,-r,U, 

Pj 

3 



(37) 



where the sum is over all the particles that lie inside radius h. 
In the GADGET code, the kernel has the following standard 
spline form: 

\2 , ^ / 



1-6U) +6(^) 
2(1-^)' 







for i < ^ < 1 



for 1< r 



(38) 



An important property of the kernel interpolant is that it 
can also be used to obtain a derivative of the reconstructed 
function, which can be simply approximated by 



, hi 



= Ej^Mrj)V^Wi\v,-v.\,h.). 
Starting from a density estimate in the form 



pi 



3 



mjW{\Yj 



(39) 



(40) 



this allows one to calculate pressure gradients, and from this, 
equations of motion for the gas elements which represent the 
Euler equations. The particular formulation for the equation 
of motion we use her e is based on the ^entropy-fo rmulation' 
of SPH discussed bv lSpringel k HernauistI (|20Q2i ). 



4.2 Obtaining anisotropic second order 
derivatives with a kernel interpolant 

Discretization of the diffusion term in the RT transfer equa- 
tion in SPH poses some difficulties. We are basically con- 
fronted with the task to find an efficient and accurate ap- 
proximation to terms of the form 



(41) 



dxsdxk ' 

where Qai3 is the product of the local Eddington tensor 
h and the relative photon density . Simply differenti- 
ating a kernel interpolant twice is not a good solution, as 
this becomes very noisy because the kernel-interpolant of 
SPH is only second-order accurate. On the other hand, the 
discretization of the Laplacian discussed by ljubelgas et al.l 
(| 20041 ) does not work either, as it only works for the isotropic 
case. 

We now describe the solution we have found for this 
problem, which basically consists of the task to approximate 
the second order partial derivatives of an element Qq;/3(x) of 
the matrix Q(x) with a kernel-interpolant. We consider a 
Taylor-series for Qc/3(xj) in the proximity of Qapi^i), i-e. 



Qa/3(Xj) - Qc./3(Xi) = VQc^f^ 



IE 



s,k dxgdx^ 



(Xj X-i)s(Xj x^)^-!- 



(42) 



O(x,-x0^ 

Let us use the short-hand notation x^j = Xj — x^ and 
Wij — W(\yLj — Xi|), where W{r) is the SPH smoothing ker- 
nel. Neglecting higher order terms, we multiply the above 
expansion with 

|x,,P ' ^^^^ 

and integrate over all Xj. Here (x^j); is the l-ih component 

>artic 

In particular, this 



with respect to the m-component of x^ 
means we have 



We now find that 



/■ 



7 

,m 



' |x,,-| 



dx^ 0, 



(44) 



(45) 



for all combinations of /c, / and m. This is because there is 
always at least one single component of x^j left so that the 
integral vanishes by symmetry. As a result, the first order 
term of our integrated Taylor expansion drops out. 

We now consider the second order term, where we en- 
counter the expression 

(Xij)s i^ij)k (^ij)l (Xij)m VF^(|Xi^ " 



■ dx, 



(46) 



There are a number of different cases. If / and m are equal, 
then s and k must also be equal, otherwise the integral van- 
ishes. So here we would have three possible contributions 
to a s,/c-sum, corresponding to the three coordinates that 
s and k can assume. If / and m are unequal, then we must 
either have s — I and /c = m, or have s — m and k — s. So 
here there are two contributions to a s,A:-sum in this case. 
Evaluating the integral Tskim for these cases gives: 



if i 



: m and s — k 



Tsklm — ^ 



I if / = m and s — k^ but s / /, 



if / / m, and s = I, k = m 
or s = m, /c = /, 



(47) 



in all other cases. 

Note that we can pick / and m freely when we multiply the 
Taylor expansion with the term (|43|) and integrate over it. 
In particular, we can also use several different choices one 
after the other and then form a linear combination of the 
results. This can in fact be used to isolate any of the second 
derivatives of the Hessian matrix of Qaf3- Let us assume for 
example that we want to calculate the second derivative of 
Qaf3 with respect to xq. Choosing / = m = 0, then the three 
choices /c = s = 0, /c = s = l and /c = s = 2 all give terms 
that contribute to the integral over the expansion. These 
are: 
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3 d^Q I 1 d^Q I 1 a^Q 

5 dx'i ^ 5 ~^ 5 



(48) 



Here Q is to be understood as Q = for brevity. Based 
on this, we can now isolate the desired partial derivative by 
forming a linear combination: 

Q(x.)-Q(x,) 



2/- 



(49) 



[2(x,^)oT^i^-,o - ^(x^^)iT^i^-,i - ^{:>c^j)2WiJ,2] dx^-. 

In a similar fashion, we can obtain a mixed partial derivative 
in the following way: 

Q(x,)-Q(x^) 



2/- 



[|(X^,■)0V^^,■,1 + |(x^,■)lV^^J,o] dx," . 



(50) 



Formulae for all other second-order partial derivatives 
can be obtained from these expressions by cyclic permuta- 
tion. Also, they are valid for each of the matrix elements 

Qaf3- 

Using these results, we can now turn to obtaining an 
expression for the sum of the second derivatives, as needed 
in the anisotropic diffusion equation. Based on the above, 
we can write the desired expression in the compact form: 



[Q(x,) - Q(xi)] ViWij 



(51) 



Here we defined a new matrix Q through the matrix ele- 
ments of the original matrix Q = (Qap), in the following 
way: 



Q = f Q - ^Tr(Q)I. 



(52) 



Inspection of this result highlights one interesting issue 
that could potentially become a numerical stability prob- 
lem in certain situations. The matrix Q is not guaranteed 
to correspond to a positive definite quadratic form when it 
is used in the SPH discretization form of equation (|5T]) . If 
the radiation transfer is very anisotropic, the matrix Q can 
contain negative diagonal elements and thus gives rise to an 
'anti-diffusive' behaviour in the discretized radiation trans- 
fer equation, where radiation is transported from a particle 
of lower radiation intensity to one with higher radiation in- 
tensity. It is not clear right away whether this will lead to 
numerical stability problems of the radiative diffusion treat- 
ment, but it could. 

In case this is a problem, one way to avoid it would 
be to somehow suppress transport of radiation opposite to 
the direction of the gradient of the radiation intensity be- 
tween a particle pair. Another way is to add in an isotropic 
component to Q such that 



(53) 



Here the idea is to make Q* slightly more isotropic, such 
that Q* becomes positive definite again. In order to guar- 
antee this, we need to assign a = |. Form equations (|52|) 
and ()53p we can see that this 'anisotropy- limited' matrix is 
then actually Q* = Q. One interpretation of this result is 
that the unmodified matrix Q mediates diffusion which is 



a mix of 2/5 of the 'correct' anisotropic diffusion and 3/5 
of isotropic diffusion. In some of our tests we will compare 
results from both formulations of the matrix. We will refer 
to Q* as the 'anisotropy-limited' tensor and to Q as the 
'fully-anisotropic' tensor. 

Note that in the case where Q is diagonal and pro- 
portional to the identity matrix, equation (|5T]) reduces 
to the isotrop i c resu lt for a scalar function derived by 
ljubelgas et al ] (|20Q4l ) for the thermal conduction problem. 
The important point about equation (|5T]) is that it involves 
only a first order derivative of the kernel function. As a 
result, it can be discretized straightforwardly in the usual 
SPH way, where the integration is replaced by a sum over 
all neighboring SPH particles within the kernel volume. 



4.3 Discretization of the anisotropic diffusion 
term in SPH 

As we have shown in detail in the previous section, kernel- 
interpolated second-order derivatives of some tensor Q^p 
can be obtained as in equation (|5T]) for a suitably defined 
modified tensor Q. Furthermore, we note the identity 



dx ys dy J 2 y dxdy s dxdy s s dxdy J 
and thus inserting equation (|51|) we obtain 



(54) 



d 



2 J ^ ' i_, r,. dx," 



(55) 



For our application, let us denote the correspondingly mod- 
ified Eddington tensor as h^-^ . We can then write down an 
SPH discretization of the diffusion part of the radiation 
transfer equation. This can be expressed as: 



'7 



dh. 
~dt 

Here 



^ ^ C Xg- [n^h, - hl^h^] V^W^J rUj 

^-^ K.^n x^. n^n ' 



1 1 
— + — 



(57) 



is a symmetric average of the absorption coefficients of the 
two particles i and j, and pij is a symmetrized density. It 
is however important to be careful about how exactly the 
symmetrizations are done in practice, because this can af- 
fect the performance of the scheme if there are particles of 
varying mass. In particular, we would like to use a formu- 
lation where the conservation of the number of photons is 
guaranteed in this case as well. If possible, we would also like 
to obtain a formulation where the effective coupling matrix 
is symmetric, because this is a prerequisite for using certain, 
particularly efficient solution methods from linear algebra, 
such as the conjugate gradient (CG) method. 

The photon conservation property is best analyzed by 
switching to variables that directly encode the photon num- 
ber of each particle. Let us define for this purpose the quan- 
tity 



(58) 
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for each particle. The real photon number of a particle is 
actually = n^m/p = mhjXu/rrip, but we here drop 
the constant Xn/rrip in the definition (|58)) . for notational 
simplicity. Multiplying Eqn. (|56|) through with rrii gives now 



at ^-^ i^ijPij x; 

Note that we can also write this as 

^N^ V--, 



where 



c :siijmihjViWij 

2 ■ 



(59) 



(60) 



(61) 



From the formulation in equation (|60p we easily see that the 
total photon number, ^.Ni, is conserved, but in general 
the matrix Wij is not symmetric. Even though other linear 
solvers may work, this would prevent us from safely apply- 
ing the CG scheme to calculate a solution for an backwards- 
Euler timestep of equation (|6Q)) , as the required implicit so- 
lution involves the inversion of a matrix that linearly de- 
pends on Wij (see below). 

To fix this problem, we also symmetrize the mass- 
weighted Eddington tensor in (|6Qp. which results in the fol- 
lowing final form of the anisotropic diffusion equation that 
we use for our numerical implementation: 

m_ 

dt 



(AT 



(62) 



(63) 



(64) 



where Wij is now redefined in a symmetric form: 
2cmij yiJjh.ijV iWij 

"Wij 2 • 

We may also include the sink term, which yields 
^ ^^Wij{Nj - Ni) - ckiNi. 

3 

This equation is still symmetric and can thus also be treated 
with the CG method, as we explain in more detail in the next 
subsection. 



4.4 Time integration of the radiative transfer 
equation 

As discussed earlier, we use an operator-split approach for 
the time integration of the radiation transfer equation, in 
which we effectively treat the time integration of the source 
term and that of the transport through anisotropic diffusion 
and the absorption through the sink terms as separate prob- 
lems. In fact, we extend the operator-split idea also to the 
hydrodynamical evolution of the system, i.e. we alternate 
the timestepping of the diffusion equation with that of the 
Euler equations that describe the dynamical evolution of the 
gas. In the following, we first discuss the time integration of 
the diffusion and sink part, which is the most complicated 
part in our scheme. 

It is well known that explicit time integration schemes of 
the diffusion equation becomes easily numerically unstable, 
unless a very small timestep is used. To ensure numerical 



stability, we therefore adopt an implicit method, namely the 
simple 'backwards Euler' scheme, which provides sufficient 
accuracy for the diffusion problem. To advance equation (|64p 
for one timestep, we therefore want to solve the equation 



(65) 



-AtcaonmK^^ , 

where A^J^^^ are the new photon numbers at the end of the 
timestep, and NJ^ are the ones at the beginning of step n. 
The last term in this equation encodes the photon loss term, 
which we also integrate implicitly. We note that the source 
term, on the other hand, can be simply advanced with an 
explicit Euler step. 

The equations (|65]) are in fact a large, sparsely popu- 
lated linear system of equations that can be written in the 
generic form 



Ax = b. 



(66) 



where A is a coefficient matrix, b is a vector of known values 
and X is a vector of unknown values. For our application, 
the components of vector b are bi = N^^ + Atsirrii, and the 
matrix elements are given by 



^ + ^AtWik 



+ Atcaonm 



Atwi 



(67) 



where the indexes i and j run over all the SPH particles. 
The solution vector of the linear problem defines the new 
photon numbers at the end of the step, Xi = N^^^ . 

There are many different approaches for solving linear 
systems of equations, but the huge size of the matrix A in 
our problem (which is equal to the particle number squared) 
makes many standard approaches that rely on storing the 
whole matrix A impractical. Fortunately, the matrix A is 
only sparsely populated because in each row only approx- 
imately ^ Nsph elements are non-zero, those that describe 
the coupling of a particle to its neighbors. Sparse systems 
of this type can often be solved well with iterative schemes. 
We use such an iterative scheme in our work, the conjugate 
gradient (CG) method. 

The conjugate gradient approach applies successive cor- 
rections to a trial solution that is used as a starting point. 
With every iteration, the solution becomes better. Since the 
corrections added in each of the steps are all orthogonal to 
each other, the rate of convergence of this method is often 
quite high, this is why we think it is a promising iteration 
scheme for the problem at hand. For reference, a deriva- 
tion of the well-known CG method is given in Appendix [B] 
and Appendix [C] However, a prerequisite for the applica- 
bility of the conjugate gradient method is that the matrix 
A is positive definite and symmetric. The symmetry is evi- 
dent from our formulation and the matrix is positive definite 
since Wij ^ 0, and thus XiAijXj ^ Vx^. 

To find a solution for the new photon number field, 
we iterate with the CG scheme until the difference between 
two successive approximations to x has dropped to a small 
percentage of |x|. Note that the expensive parts in the calcu- 
lation of one iteration are the matrix- vector multiplications. 
For each particle, they reduce to sums over all of its SPH 
neighbors, which is equivalent to an ordinary SPH loop, sim- 
ilar in computational cost to, e.g., the SPH density estima- 
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tion. Since in our parallel code some of the SPH neighbours 
of a particle can be stored on other processors, this step also 
involves communication. 

In practice, we start the iteration with the current pho- 
ton distribution for x, which is usually a fairly good starting 
point, since the expected solution for the photon distribution 
does not differ significantly from the previous state, except 
in the vicinity of the sources. In the absence of a radia- 
tion field we set the vector to zero. The closer the guessed 
values for the vector x are to the solution, the faster the 
algorithm converges. One could also start from a random 
photon distribution, which will not affect the solution, but 
will slow down the algorithm and is therefore not a desirable 
choice. The number of iterations required to reach conver- 
gence depends on the condition number of the matrix A, 
where the condition number is defined as A max/Amin, the ra- 
tio of the largest to the smallest eigenvalue of the matrix. 
A large condition number slows down the convergence rate 
of iterative solvers of linear systems of equations. However, 
the number of required iterations can be reduced by pre- 
conditioning the matrix A. For this purpose we employ the 
simple Jacobi preconditioner. More specifically, we modify 
our matrix equation as follows, 

C"^Ax = C"^b, (68) 

where the matrix C is the Jacobi preconditioning matrix. It 
is defined as 



(69) 



(70) 



Cij — AiiSij ^ 

and has the rather simple inverse form 

1 _ 6ij 

Applying the Jacobi preconditioner to the matrix A basi- 
cally means to divide the matrix A by its diagonal, which is 
simple to implement and to parallelize. While the associated 
reduction of the condition number improves the convergence 
speed, it would be desirable to find still better precondition- 
ers that are more effective in this respect. 



4.5 Implicit time integration of the chemical 
network 

The chemical network we need to solve is described by equa- 
tion (|35p . Even though it contains only a simple time deriva- 
tive, the numerical integration can be rather tricky because 
very large changes in the number of ionising photons can 
suddenly occur from timestep to timestep. A stable integra- 
tion with a simple explicit scheme therefore requires a very 
small timestep, of the other of 10~^ the typical dynamical 
timestep of the simulation. This would make the solution of 
the chemical network very expensive, and it is highly desir- 
able to find another method that speeds up the computation. 

To this end we use an implicit scheme for the evolu- 
tion of the chemical network. Our variant of this approach 
discretizes equation (|35)) in time as follows (here the upper 
index n denotes the integration time step): 



^HII 



At 



- ann {nun) , 



(71) 



where we have substituted he = hmi. This quadratic equa- 
tion is easy to solve. We have found this scheme to be quite 
accurate and robust even for moderately large timesteps. 



4.6 Photoheating and cooling with a 
multi-frequency scheme 

Since our code tracks monochromatic radiation, we need to 
provide a multi- frequency description that allows photoheat- 
ing of the irradiated gas. Such approaches have already been 
introduced bv lMihalas k Weibel MihalasI (1 19841) and a re of- 
ten used in the literature (e.g. Aubert &: Tevssierll2008l V The 
frequency spectrum of the source is approximated by a black 
body spectrum, which allows the use of a frequency indepen- 
dent photoionisation cross-section a defined as 



Jo 



diy — ; X 



dv 



4:7VJu 



(72) 



Similarly, a frequency averaged photon energy e (above the 
hydrogen ionisation energy) can be defined. Thus the hy- 
drogen photoheating rate from equation ([8]) can be approx- 
imated by 



(73) 



where c is the speed of light and is the number density 
of ionising photons. 

In order to evolve the temperature of the gas correctly 
we consider cooling processes such as recombination cooling, 
collisional ionisation cooling, collisional excitation cooling 
and Bremsstrahlung cooling. All rates have been adopted 
from Cen (1992) and are summarized in Appendix [Al As a 
combination of heating and cooling terms the temperature 
of each particle is then evolved as follows 



+ At- 



[r-A], 



(74) 



where the upper index n denotes the timestep of the sim- 
ulation, /cb is the Boltzmann constant, nu is the hydrogen 
number density, F is the heating rate and A is total cooling 
rate, evaluated for the new temperature at the end of the 
step. 

4.7 The source terms 

As given in equation (|34p . the source term in the RT equa- 
tion is in our formulation. This represents the number of 
photons emitted per unit time and per hydrogen atom. In 
our cosmological applications we usually want to represent 
this term by the stars that have formed in the simulation. 
To more accurately account for the short-lived massive stars, 
we can also use star-forming SPH particles as sources, con- 
verting their instantaneous star formation rates into an ion- 
ising luminosity. Alternatively, we can also consider harder 
sources like Active Galactic Nuclei (AGN), if they are fol- 
lowed in the simulation. The time integration of the source 
function is unproblematic, and can be done with a simple 
explicit scheme on the dynamical timestep of the simulation. 
We only need to assign the source luminosities in a conser- 
vative fashion to the nearest SPH particles, provided they 
are not already gas particles anyway. 

In some of our test problems considered in Section O the 
source term is prescribed as a certain number of ionising 
photons per second, independent of the mass of the source. 
In this case we have 



nu V 



N ^ 



(75) 
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Figure 1. Eigenvectors of the Eddington tensor for a single source (left panel), and for two sources (right panel), calculated with our 
tree-code extensions of the GADGET-3 code. Both vector fields match the expectations. Note that the directions of the vectors can be 
turned 180° without affecting the direction of the transport of radiation, an effect due to the symmetric nature of the Eddington tensor. 



where 1/ = m/p is a measure of the volume of the SPH 
particle with mass m that hosts the source. If we choose to 
distribute the source function to more than one particle, we 
use the 'scatter' approach. In this case the source function 
has the form 

TTlj 



(76) 



where the sum is over the SPH neighbors of the emitting 
particle. 



4.8 Eddington tensor calculation 

An important quantity in our formulation of the RT problem 
are the local Eddington tensors, which we estimate based 
on an optically thin approximation, as defined by equations 
(pO)) and ([2T]) . The dependence of the contribution of 
each source suggests a calculation method similar to that 
of gravity - via a hierarchical tree algorithm. To this end 
we extend the gravitational tree code in GADGET-3 with 
additional data structures. For each node of the tree, we also 
calculate the total luminosity and the luminosity-weighted 
centre-of-mass. Depending on whether we consider the star 
particles as sources, the star- format ion rate of gas particles, 
or also black hole particles, this can involve different types 
of particles. 

The individual elements of the Eddington tensor are 
then computed by walking the tree in a way exactly analo- 
gous to the procedure for calculating gravitational forces. If 
a tree node appears under a small enough angle as seen from 
the point of interest, all its sources can be represented by 
the total luminosity of the node. Otherwise, the tree node is 
opened, and the daughter nodes are considered in turn. As 
a result, we obtain a multipole approximation to the local 
radiation pressure tensor, with a typical accuracy of ^ 1%, 
provided a similar node-opening threshold value is used as 
for collisionless gravitational dynamics. However, since for 
the Eddington tensor only the direction of the radiation 



pressure tensor ultimately matters, the final accuracy of the 
Eddington tensor is even better, and more than sufficient 
for our purposes. A very important property of this calcu- 
lational method is that its overall speed is essentially inde- 
pendent of the number of sources that are present since the 
calculation is done simultaneously with the tree-walk that 
computes gravity, and involves only few additional float- 
ing point operations. This is quite different from the widely 
employed ray-tracing or Monte- Carlo schemes for radiation 
transfer, where the calculation cost may scale linearly with 
the number of sources. 

As an example. Figure [T] shows the eigenvectors of the 
Eddington tensor for two different source configurations, cal- 
culated with our modified version of the GADGET-3 code. 
The left panel is for a single source. The vectors point radi- 
ally outward or inward from the source, as expected. Note 
that the directions of the vectors can be turned by 180° with- 
out changing the radiation transport, because the tensor is 
symmetric. The panel on the right hand side shows the vec- 
tor field around two sources, in the plane of the stars. The 
field is a dipole in this case and matches the expectations 
well. 



4.9 Flux limited diffusion 

If is the density of photons, then the maximum photon 
flux / that can occur is limited by the speed of light to 
/ = cn^. This physical limit for the possible photon flux 
can in principle be violated under certain conditions when 
the diffusive approximation to the photon flux. 



1 Idn^h'' 



(77) 



is used. In treatments of radiative transfer in the isotropic 
diffusion approximation one therefore often invokes so-called 
flux limiters that are designed to enforce the condition 



(78) 
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Figure 2. Resolution comparison of spherically averaged ionised and neutral fractions as a function of radial distance from the source, 
normalized by the Stromgren radius tq. Results are shown for an anisotropy-limited (left panel) and fully anisotropic (right panel) 
Eddington tensor formulation at time t = 500 Myr ~ 4trec- The compared resolutions are 8^, 16^, 32^, 64^ and 128^ particles. The black 
lines show the analytical solution, integrated radially outward from the source as in equation (|87|) . In both formulations the accuracy 
clearly increases with resolution and saturates when the highest number of particles is reached. In all cases, the final Stromgren radius 
agrees with the analytical predictions. The anisotropy-limited Eddington tensor formalism gives very accurate predictions for the ionised 
fraction in the regions outside the Stromgren radius (r > rs), but fails to give a correct value for the inner parts of the ionised regions. 
The fully anisotropic Eddington tensor formalism, however, predicts accurate values in both regions. 




Figure 3. Resolution comparison of the scatter (grey areas) of the spherically averaged ionised and neutral fraction as function of radial 
distance from the source, normalized by the Stromgren radius rs- Results are shown for an anisotropy-limited (left panel) and fully 
anisotropic (right panel) Eddington tensor formulation at time t = 500 Myr ~ 4frec- The compared resolutions are 16^, 32^, 64^ and 
128^ particles. The black lines show the analytical solution, integrated radially outward from the source as in equation (|87|) . The orange 
lines show the spherically averaged neutral fraction and the violet lines the spherically averaged ionised fraction. The range of the scatter 
does not change significantly with resolution, since it is due to the intrinsic diffusive nature of SPH and the inaccuracies of the SPH 
density estimate. The highest resolution simulation has no scatter due to the high accuracy of the density estimates. The scatter in the 
fully anisotropic Eddington tensor formalism is larger due to the larger anisotropy in the diffusion terms of the RT equation. 



by damping the estimated flux when needed. ^ _ I 1 ^^g^^ 

In our anisotropic diffusion treatment, we have also im- 

plemented an optional limiter that serves the same purpose. 

1 ,1 . n . • . u u 1 f We then define a flux limiter of the form 

We observe the maximum flux constramt with the help of 

a parameter which is a function of the gradient of the i + 0.1i? 

photon density ^(^) = TToI^ToOLR^ ' ^^^^ 
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Figure 4. I-front expansion for an anisotropy-limited (left panel) and a fully anisotropic (right panel) Eddington tensor formulation 
with 16*^, 32"^, 64*^ and 128"^ particles. The dashed line is the exact solution obtained from equation (|87|) and the solid one as obtained 
from equation (j83|) . Both results agree very well with the analytical predictions, but the fully anisotropic Eddington tensor formulation 
shows better results and smaller relative error for the I-front position. 



where X{R) ^ as ^ oo. The detailed 
lytic expression used for the flux limiter is 
as it ensures a smooth transition between 
states. We have chosen this form since it 
other nu merical RT codes, fo r exam ple, a 
used by IWhitehouse k Bate! (|2004l ). The 
then introduced into the diffusion part of 
follows: 



^ d 



K, dxi 



form of the ana- 
arbitrary as long 
the two limiting 
is widely used in 
similar version is 
flux limiter A is 
equation ()34p as 



(81) 



However, we note that super luminal propagation of 
photons is usually not a problem in the ionisation prob- 
lems we are interested in. Here the speed of the ionisation 
fronts is not limited by the speed of light, but rather by the 
luminosity of the sources and the density of the absorbing 
medium. Nevertheless, the flux limiter can also be usefully 
employed as a means to control the behavior of the ionisa- 
tion front (I-front) propagation in dense media. Due to the 
specific dependence of the R parameter, the photon prop- 
agation can effectively be limited in high-density regions, 
where the intensity gradient becomes very large. 



4.10 Notes on the performance of the code 

An important consideration in the development of RT codes 
is the calculational cost of the implemented scheme, as this 
determines whether the method is sufficiently fast to allow 
a coupling with hydrodynamic simulation codes. It is diffi- 
cult to make general statements about the computational 
cost of our new RT scheme as this depends strongly on 
the particular physical problem that it is applied to. Ar- 
guably of most interest is a comparison of the speed of the 
method to other parts of the simulation code when applied 
to the problem of cosmological structure formation with a 
self- consistent treatment of cosmic reionisation. We here re- 
port approximate numbers for the speed of our new method 



in this situation, based on work in preparation that studies 
this problem. 

First we note that in our implementation, when gravity 
is also integrated, the computation of the Eddington ten- 
sor incurs negligible costs, since it is done together with the 
gravity. Therefore, we find that in our implementation the 
increase of the computational cost with respect to the other 
relevant code modules (gravity, SPH density and SPH hy- 
drodynamical forces) is primarily a function of the number of 
iterations required at each timestep to construct the implicit 
solution of the anisotropic diffusion equation. One iteration 
is approximately as costly as one SPH hydro computation, 
and the average number of iterations required ranges from 
typically 10 up to 200 in the most extreme cases. However, 
the RT equation is integrated only on a relatively coarse 
timestep, whereas the hydrodynamics is done also for many 
more smaller sub-steps. This reduces the effective cost of 
the whole RT calculation to several times the total cost of 
the hydrodynamics calculations. In practice, we measure a 
slow-down of the simulation code by a factor of order of 2 
to 5 when the radiative transfer is included. While this is 
non-negligible, it does not seriously impact the ability to 
carry out large cosmological simulations. We also note that 
further optimizations in the radiative transfer algorithms, 
perhaps through an improved pre-conditioner, may reduce 
the cost of the RT calculation in the future. 



5 TESTING THE CODE 

In the following we present several basic tests for our new 
radiative transfer code. Where possible, we compare our re- 
sults with analytical solutions or with results from other sim- 
ulatio ns (from the RT code comparison study bv llliev et al.l 
l2006l ). In section O we discuss our results for the classic 
test of the isothermal expansion of an ionised sphere in a 
homogeneous and static density field. In section 15.1.11 we 
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study the effects of different timesteps on the accuracy of 
our simulations, while in section [5. 1.21 we evolve two nearby 
sources with interacting ionised spheres. Then in section [5121 
we repeat the single ionised sphere expansion test, but this 
time allowing the temperature to evolve. In section 15.31 we 
present a shadowing test, where a dense clump is placed in 
the way of a plane-parallel I-front. Finally, in section 15.41 
we evolve the radiation transport in a static cosmological 
density field. 



5.1 Isothermal ionised sphere expansion 

The expansion of an I-front in a static, homogeneous and 
isothermal gas is the only problem in radiation hydrody- 
namics that has a known analytical solution and is therefore 
the most widely used test for RT codes. A monochromatic 
source emits steadily photons with energy hu = 13.6 eV 
per second into an initially neutral medium with constant 
gas density nn- Then the Stromgren radius, at which the 
ionised sphere around the source has reached its maximum 
radius, is defined as 



rs 



3N^ 



1/3 



(82) 



where as is the recombination coefficient. This radius is 
obtained by balancing the number of emitted photons by 
the number of photons lost due to recombinations along a 
given line of sight. If we assume that the I-front is infinitely 
thin, i.e. there is a discontinuity in the ionisation fraction, 
then the expansion of the Stromgren radius can be solved 
analytically and the I-front radius n is given by 



n = rs[l - exp(-t/trec)]^'^^, 



where 



(83) 



(84) 



is the recombination time and is the recombination co- 
efficient. 

The neutral and ionised fraction as a function of radius 
of th e stabl e Stromgren sphere can be calculated analytically 
(e.g. lOsterbrock &; Ferlandll2006l ) from the equation 
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where riHi is the neutral fraction, riHii is the ionised fraction 
and 
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Moreover, considering spherical symmetry and a point 
source we can solve analytically for the photon density radial 
profile n^(r), yielding 
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(87) 



From this we obtain ionised fraction profiles nHii(^) for the 
whole evolution time. 

The Stromgren radius obtained by direct integration of 
equation ()85p differs from the one obtained from equation 
(|82|) because it does not approximate the ionised region as a 



sphere with constant, but with varying ionised fraction. We 
compare our results with both analytical solutions. 

For definiteness, we follow the expansion of the ionised 
sphere around a source that emits Nj = 5 x 10^^ photons s~^ . 
The surrounding hydrogen number density is nu = 
10~^ cm~^ at a temperature of T = 10^ K. At this temper- 
ature, the case B recombination coefficient is as = 2.59 x 
lO""*^^ cm^ s""*^. Given these parameters, the recombination 
time is tree = 125.127 Myr and the expected Stromgren ra- 
dius is rs = 5.38 kpc. We impose periodic boundary condi- 
tions in order to make sure that the density field is effectively 
infinite and uniform. We note that this does not affect our 
RT calculation since the Eddington tensor is computed non- 
periodically. 

We present results from the fully-anisotropic and the 
anisotropy-limited Eddington tensor formalism simulations 
by comparing first the final state of the ionised sphere 
and then the evolution of the I-front. Figure [2] shows the 
spherically averaged ionised and neutral fraction as a func- 
tion of radial distance from the source, normalized by the 
Stromgren radius rs, at time t = 500 Myr ^ 4 tree- In this 
case we compare the resolution effects on the accuracy of 
our numerical predictions by using simulations with 8^, 16^, 
32^, 64^ and 128*^ particles, corresponding to mean spa- 
tial resolutions of 2.0 kpc, 1.0 kpc, 0.5 kpc and 0.25 kpc. 
In both formalisms the accuracy increases with resolution 
and the profiles converge for the higher resolutions. We also 
note that the anisotropy-limited Eddington tensor formal- 
ism gives very accurate predictions for the ionised fraction 
in the regions outside the Stromgren radius (r > rs), but 
fails to give correct values in the inner parts of the ionised 
sphere. The fully anisotropic Eddington tensor formalism, 
on the other hand, predicts the ionisation state in both re- 
gions quite accurately. 

We compare also the scatter of the ionised and neutral 
fraction profiles. Figure [3] shows the scatter (gray areas) of 
the spherically averaged ionised and neutral fraction profiles 
for four different resolutions (16^, 32^, 64^ and 128^ parti- 
cles) at the end of the I-front expansion. All results agree 
well with the analytical radius of the Stromgren sphere. The 
range of the scatter does not change significantly with res- 
olution, since it is due to the intrinsic diffusive nature of 
SPH and the inaccuracies of the SPH density estimate. This 
means that we obtain a density scatter of about 0.01% and 
thus introduce fluctuations in the gas density, which result 
in fluctuations in the hydrogen densities of the SPH particles 
and thus of the ionised and neutral fractions. The SPH den- 
sity scatter in the 128^ particle simulation is zero (thanks 
to the use of a Cartesian grid - but note that in real-world 
dynamical applications some density scatter is unavoidable) 
and thus there is no scatter in the ionised and neutral frac- 
tions. Moreover, the scatter in the fully anisotropic Edding- 
ton tensor formalism simulations is larger than in the other 
formalism simulations due to the larger retained anisotropy 
in the diffusion term of the RT equation. 

The evolution of the I-front expansion is shown in Fig- 
ure 2] comparing the two formalisms at different resolutions. 
The results from both formalisms agree very well with the 
analytical predictions and the accuracy increases with res- 
olution. The fully anisotropic Eddington tensor formalism 
simulations show better results and smaller relative error 
for the I-front position. In both formalisms the error stays 
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Figure 5. Slices of the neutral fraction through the ionised sphere at the position of the source. The upper row shows the results from the 
anisotropy-limited Eddington tensor formulation and the lower row from the fully anisotropic Eddington tensor formalism. Five different 
resolutions are compared: 8^, 16^, 32^, 64^ and 128^ particles. The contours mark neutral fractions of tt-hi = O-Q? 0-5, 0.1, 0.01 and 0.001. 
The white circles give the radius of the Stromgren sphere. The geometrical distribution of the SPH particles affects the shape of the 
ionised sphere. The spheres are elongated in the x- and y-direction of the Cartesian grid, where the particle spacing is smaller and the 
SPH-kernel interpolant leads to slightly different couplings as in diagonal directions. This effect is stronger for the anisotropic Eddington 
tensor formulation. 



within 5% of the analytical solution and traces the analytical 
result obtained by direct integration of equation (|87)) . 

However, we note that the geometrical distribution of 
the SPH particles we used in our simulations introduces 
slight deviations from perfect sphericity into the shape of 
the ionised region. This reflects the Cartesian grid of parti- 
cles used for these tests, an effect that can be clearly seen 
in the shapes of the ionised regions displayed in Figure O 
The spheres are elongated in the x- and ^-directions of the 
Cartesian grid, where the particle spacing is smaller and 
the SPH-kernel interpolant weights the nearest neighbours 
slightly differently than in off- axis directions. This discrete- 
ness effect is stronger for the anisotropic Eddington tensor 
formulation. 

Considering all our results that compare the limited and 
fully anisotropic Eddington tensor formalisms, we use in all 
our further tests and simulations only the fully anisotropic 
formulation because it shows more accurate results. It turns 
out that this formulation is also robust, i.e. it does not show 
stability problems due to its 'anti-diffusive' terms when used 
in conjunction with an implicit solver, at least we have not 
experienced such problems in our test calculations. 

5.1.1 Timestep comparison 

In order to test the accuracy of our RT scheme we perform 
simulations of 64"^ particle resolution with different fixed 
timesteps: At = 0.05, 0.5, 5and50Myr. Applying the von 
Neumann stability criterion for an explicit integration of the 
diffusion part of our RT equation (|19|) , we find a bound of 
the timestep equal to 



where Ax is the mean spatial resolution, c is the speed of 
light and k = nm cf is the absorption coefficient. At the I- 
front the assumed neutral fraction is nni =0.5 and thus the 
absorption coefficient is />: = 3.15 x 10~^^ cm~^, resulting in 
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Figure 6. I-front expansion as a function of the recombination 
time. Results are shown for four different simulation timesteps: 
At = 0.05, 0.5, 5and50Myr. The dashed line is the exact solu- 
tion obtained from equation (|87|) and the solid one is obtained 
from equation (183 1) . The smallest timestep simulation agrees very 
well with the analytical solution. As the timestep increases, the 
results in the early phases of the expansion become more inaccu- 
rate. However, after about two recombination times, the I-front 
radius catches up with the analytical solution. The simulation 
with timestep At = 50Myr is very inaccurate, but in the end of 
the expansion the I-front radius is still within 5% of the analytical 
solution. 

an upper limit for the timestep At ^ 10""^ Myr. However, 
this limit on the timestep is only a reference point for our re- 
sults. Because we use an implicit scheme that is stable for all 
timestep sizes, we are fortunately not bound by this timestep 
limit and can in principle use much larger timesteps, subject 
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Figure 7. Ionised fraction in a plane of two equally luminous 
sources. The positions of the sources (8kpc apart) are marked 
with black crosses. The snapshot is taken at time t = 500 Myr, 
when the expansion of the ionised regions has stopped, and the re- 
gion where the Stromgren spheres overlap is approximately 3 kpc 
wide. There is a clear elongation along the axis c onnec ting the 
two sources, as also described bv lGnedin Sz Abe] (|20Qll V Figure 
[8] shows the evolution of the I- front along the aligned (red) and 
perpendicular (blue) directions with respect to the axis through 
the sources. 



only to the condition that the final accuracy reached is still 
acceptable. 

We compare results for our simulations with the four 
different timestep sizes in Figure (6] The smallest timestep 
that we use is ten times as big as the analytical upper limit 
for an explicit scheme, yet its results agree perfectly with the 
analytical solution. For the other simulations the timestep 
sizes increase progressively by factors of 10, and the numer- 
ical results start to deviate from the analytical calculation. 
The changes are largest in the early phase of the I-front ex- 
pansion. As the source "suddenly" switches on, very small 
time steps are needed in order to achieve good accuracy in 
the beginning, where the gradients in the photon density are 
very large. But later, after a couple of recombination times, 
the numerical results approach the analytical solution even 
for coarse timesteps and follow it until the expansion of the 
I-front ends. The At = 50 Myr simulation is initially very 
inaccurate, but note that its Stromgren radius is still within 
5% of the analytical result. Therefore, our method manages 
to essentially produce correct Stromgren radii of the ionised 
spheres for all considered timesteps. 



5.1.2 Two nearby sources 

In our next test we follow the expansion of ionised regions 
around two nearby sources, where we expect to see inac- 
curacies due to the optically thin assumption used for es- 
timating the Eddington tensors. Both sources emit = 
5 X 10^^ photons s~^ and are 8 kpc away from each other. The 
number density of the surrounding static and uniform hydro- 
gen gas is riH = 10~^ cm~ ^, at a temperature of T = 10^ K. 
From tests conducted bv iGnedin Abel (|200lh we expect 
that the ionised regions are not spherical, but rather elon- 
gated along the axis through the sources. This effect results 
from the calculation of the Eddington tensor, whose values 
along the symmetry axis are estimated high and boost the 
diffusion in this direction, while reducing it in the perpen- 
dicular direction. 

In Figure [3 we show the neutral fraction for this test in 
a slice in the plane of the sources, taken at time t = 500 Myr 
when the expansion of the regions has stopped. The expected 
elongated shape of the ionised regions is clearly visible. In 



2.0 



1.8 




0,8 



0.6 



sphericolly overoged 
. perpendiculor 
. oligned 



100 



200 300 
t [Myr] 



400 



500 



Figure 8. I-front expansion around a source in a double sys- 
tem, normalized by the analytical solution from equation (|83|) . 
as a function of time. The red line is the I-front in the direction 
aligned with the two sources and the blue line is for the orthog- 
onal direction. The green lin e shows the spherically averaged I- 
front position. As observed bv lGnedin Abell (|200ll ). the ionised 
spheres are elongated along the axis of the two sources and com- 
pressed in the perpendicular direction. The spherically averaged 
I-front position is within 20% of the analytical expectation. 



Figure [8] we show the time evolution of three character- 
istic radii of the expanding ionised regions: one radius is 
measured in a direction aligned with the axis through the 
sources, one is measured perpendicular to it, and the third is 
a spherically averaged radius. We note that we do not expect 
the radii to match the analytical prediction from equation 
(|83|) exactly since the approximations there are valid only 
for a single ionised region expansion, but we here use the 
obtained value to compare the expansion of the ionised re- 
gions around two nearby sources. As expected, the aligned 
radius is always larger than the analytical result, while the 
perpendicular radius is smaller. However, the spherically av- 
eraged radius of the expanding region stays within 20% of 
the analytical value. We conclude that the optically thin 
approximation to estimate Eddington tensors can in certain 
situations introduce errors in the shapes of ionised bubbles, 
but these errors should be quite moderate or negligible in sit- 
uations where the Eddington tensor is dominated by a bright 
nearby source, which is probably generic in many scenarios 
for cosmological reionisation. 



5.2 Ionised sphere expansion with varying 
temperature 

In this test we use the same setup as in section 15. 1[ but 
we initialize the gas temperature with To = 10^ K and let it 
evolve due to the coupling to the radiation field. We further- 
more approximate the recombination coefficient as with 



aB(T)=2.59xlO-"(^)" 
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cm s , 



(89) 
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Figure 9. Time evolution of temperature (left panel) and ionised and neutral fraction (right panel) profiles of an expanding ionised 
sphere for two different resolutions: 16^ (blue) and 64*^ (red) particles. Results at times t = 10, 100 and 500 Myr are shown in solid, dotted 
and dashed lines, respectively. The temperature profiles inside the ionised sphere converge for both resolutions at all times. The position 
of the I- front agrees well for both resolutions at all times. 
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Figure 10. Volume fraction of the temperature for a comparison between CRASH, OTVET and GADGET at three different times 
t = 10, 100 and 500 Myr (left to right). The results of GADGET and OTVET are comparable, but CRASH has a harder spectral 
distribution and employs multiple frequency bins and thus gives slightly different results. Differences between GADGET and OTVET 
may also be due to different resolutions (OTVET has 0.05 kpc, while for GADGET ~ 0.25 kpc). 



and assume a black body spectrum of temperature 10^ K 
for the source, setting the parameters from equation (ff3)) to 
a = 1.63x10"^^ cm^ and e = 29.65 eV. Then we evolve equa- 
tion ()74p for every particle, at every time step, considering 
photoheating, recombination cooling, collisional ionisation 
cooling, collisional excitation cooling and Bremsstrahlung 
cooling. In this way we test a realistic expansion of an ionised 
sphere around a single source. 

We test our scheme with two different resolutions: 
16^ and 64^ particles. The time evolution of the temperature 
profile in both tests is shown in the left panel of Figure (9] 
The temperature close to the source rises to approximately 
2 X 10^ K and then drops down to 10^ K outside the ionised 
region. From left to right, the profiles are shown at three 
different times: t = 10, 100, and 500 Myr. The results from 
the simulations converge inside the ionised sphere. Outside 
the ionised region the low resolution simulation produces. 



as expected, a smaller slope of the temperature drop further 
from the ionised sphere radius. In the right panel of Figure 
[9l we compare the neutral and ionised fractions at the same 
times. Both resolutions converge at the radius of the ionised 
sphere. 

In order to verify our results we co mpare them with re - 
sults obtained with the codes CRASH f' Maselli et £11120031 1 
and OTVET (Gnedin & Abel 2001), as summarized by 
llliev et al.l (|2QQ6l l in the RT Code Comparison Project. In 
Figure [lOl we show a comparison of the volume fraction of 
the temperature and in Figure [11] we present a compari- 
son of the ionised volume fractions, at three different times: 
t = 10, 100, and 500 Myr. The temperature volume fractions 
are all somewhat different due to the very different heating 
schemes employed. The results of GADGET and OTVET are 
comparable, but CRASH uses a harder spectral distribution 
and multiple frequency bins and thus gives different results. 
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Figure 11. Volume fraction of the ionised fraction in a comparison between CRASH, OTVET and GADGET for three different 
times t = 10, 100 and 500 Myr (left to right). With increasing time, GADGET produces larger strongly ionised volume and smaller 
intermediately ionised volume than the other codes, mirrored in the larger gradient of the ionised fraction that it produces. These 
deviations are likely due to different source spectra treatments, temperature structures of the ionised spheres, and different numerical 
resolutions. 



With increasing time GADGET produces a larger strongly 
ionised volume, and smaller intermediately ionised volume 
fraction than the other codes, which is also mirrored in the 
larger gradient of the ionised fraction that it produces. These 
differences are due to the different treatments of the source 
spectra that the codes employ, which are in general difficult 
to compare. Deviations might also be due to the different 
temperature structures of the ionised spheres and the dif- 
ferent resolutions of the codes (OTVET and CRASH used 
0.05 kpc, GADGET only - 0.25 kpc). 

5.3 Shadowing by a dense clump 

As a further test problem, we consider the interaction of a 
plane-parallel front of ionising photons with a uniform dense 
cylinder of neutral gas. The setup of our problem consists 
of a box with dimensions {x,y,z) = (40 kpc, 10 kpc, 20 kpc). 
One side of which (the xy-p\a,ne) is aligned with a plane of 
stars that produce the ionising photons. A dense cylinder 
of gas is located 5 kpc from this sheet-like source and has 
a radius rc = 2.5 kpc. The axis of the cylinder is oriented 
parallel to the z-axis of the box. The particle resolution is 
{Na:,Ny,Nz) = (256,64,8). The hydrogen number density 
in the cylinder is 10^ times the surrounding density of nu = 
10~^cm~^. There are 512 stars and each of them emits — 
1.2 X 10^^ photons s"\ 

We first present results obtained with an Eddington ten- 
sor that mimics a plane-parallel I-front, of the form 

~ / 1 \ 

h = , (90) 
V / 

which represents photon transport only in the x-direction. 
We also consider two other cases where the Eddington tensor 
is calculated differently. In the first case we use the OTVET 
approximation, where the tensor is computed assuming all 
gas is optically thin. In the second case we account for the 
fact that the dense clump is optically thick by enforcing that 
no radiation is transported into the shadowed region. We 
achieve this by setting the Eddington tensor in the shadowed 
area to zero, such that the product K'^ Ji, vanishes. Note that 



for a vanishing radiation pressure tensor, the trace condition 
for the Eddington tensor needs not to be fulfilled. This sec- 
ond case is only meant to produce the expected solution 
with a sharp shadow. 

Figure [12] shows the ionised fraction in a cut through the 
simulation volume, in the plane z — ^box/2, at three different 
times: t = 0.5, 1, and 2 Myr. In the case of the OTVET 
approximation (first row) , the dense cylinder fails to produce 
a sharp shadow. Instead, the radiation also diffuses around 
the obstacle, and is propagated eventually also in directions 
different from the x-direction, albeit more slowly. However, 
if we account for the fact that the dense clump is optically 
thick (second row), a clear and sharp shadow is produced, as 
expected. In Figure [131 we show the position of the I-front 
with respect to the centre of the dense clump as a function 
of time. It is clear that the I-front moves faster in the case 
where the Eddington tensor is approximated closer to the 
analytical case, i.e. it is zero in the shadow area. 

The inability to create a sharp shadow is a limitation 
of the OTVET approximation implemented in SPH, as we 
have shown above. This limitation of codes using a moment 
method together with an OTVET approx imation has al- 
ready been noted by seve r al fflo ups, e.g. iGnedin Abel 
(|200lh : lAubert k Tevssierl (|2008h . It appears that despite 
our attempts to fully account for the anisotropic diffusion, 
the diffusion operator we have derived remains quite dif- 
fusive in SPH. This is simply a consequence of the non- 
vanishing coupling between particles with separation vectors 
not perfectly aligned with the direction of radiation propa- 
gation (here the x-axis). Unfortunately, this deficiency may 
also have other detrimental effects besides just slowing down 
the expansion of the I-front itself. We acknowledge that this 
can be an important limitation of our scheme for certain ap- 
plications, especially when shielding is common and shadow- 
ing is important. Nevertheless, this limitation can influence 
the morphology of the reionisation, but properties such as 
redshift and duration of the process, as well as temperature 
evolution of the gas and the radiation field, should still be 
accurate. 
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Figure 12. Ionised fraction (yellow) in a slice through the middle of the simulation box. Time increases from left to right and the three 
different columns of snapshots are taken at t = 0.5, 1, and 2Myr. In the upper row the Eddington tensor has been approximated with 
the OTVET scheme. In the second row the Eddington tensor has been set to zero in the shadowed region, given that the dense clump is 
optically thick, which should yield the correct result. The optically thin Eddington tensor approximation fails to produce a sharp shadow 
since radiation diffuses away from the plane-parallel front direction. In the optically thick approximation, however, a sharp shadow is 
obtained, implying that the failure of our method to produce a sharp shadow can also be blamed on using a non-vanishing Eddington 
tensor in the shadowed region. 




Figure 13. Position of the I-front, relative to the centre of the 
dense clump, as a function of time. The I-front moves faster in the 
case where the Eddington tensor is computed accounting for the 
optically thick clump (blue line). In the case where the OTVET 
scheme is used (red line), the I-front expands slower. The dashed 
line shows the position at the end of the clump, where the differ- 
ence between the expansion rates begins to grow. 



5.4 Static cosmological density field 

In our final and most demanding test calculation we fol- 
low hydrogen ionisation in a realistic cosmological density 
field, which is taken to be static for simplicity. In order to 
compare our results with those o f the cosmologic al radia- 
tive transfer comparison project (jlliev et "aL 

\mm 

we use 

the same cosmological box parameters and assign sources 
in the same way. The box with size 0.5 /i~^comoving Mpc is 
evolved with a standard ACDM model with the following 
cosmological parameters: — 0.27, = 0.043, h — 0.7, 



Figure 14. Time evolution of the volume- averaged ionised frac- 
tion in the whole simulation box. in comparison to CRASH (blue) 
and FTTE (green), GADGET (orange) produces lower ionised 
fractions at earlier times and higher ionised fractions at later 
times. The overall trend in the increase of the ionised fraction 
is the same, and the speed of ionisation first accelerates and then 
decelerates. Differences between our results and the ones from 
the Comparison Project are in part also due to differences in the 
initial conditions, as we had to transform the grid-based density 
field to an SPH realization. 



until redshift z = 9. The density field at this point is con- 
sidered for our further analysis. 

The source distribution is determined by finding halos 
within the simulation box with a FOF algorithm and then 
assigning sources to the 16 most massive ones. The photon 
luminosity of the sources is 

N. = f.^^, (91) 
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where M is the total halo mass, ts = 3Myr is the lifetime 
of the source, rrip is the proton mass and fj = 250 is the 
number of emitted photons per atom during the lifetime 
of the source. We find that the total source lu minosity in 
our si mulated box agrees well with the one from llhev etajj 
For simplicity we also set the initial temperature of 
the gas to 100 K throughout the whole box. 

We found that simply mapping the grid cells onto a 
Cartesian mesh of SPH particles with different masses intro- 
duces large noise into our RT calculation, due to the large 
variations in the mass of neighbouring particles. It is there- 
fore not straightforward to translate the mesh-based data of 
the code comparison project into an equivalent SPH real- 
ization, and we therefore needed to created our own initial 
conditions. We note that simple methods to create an equal 
particle mass SPH realization from the given grid cells, e.g. 
through random sampling, tend to introduce large amounts 
of Poisson noise and wash out extrema in the density field. 

The evolution of the total volume- averaged ionised frac- 
tion in the box is very similar for GADGET and the Compar- 
ison Project codes, as shown in Figure [Ml In the beginning 
of the simulation the total ionised fraction rises rapidly and 
then the increase decelerates. GADGET produces an overall 
lower ionised fraction until approximately t = 0.2 Myr and 
higher one at later times. This mismatch is in part certainly 
caused by the morphological differences in the initial condi- 
tions that introduce different clumping properties of matter 
and therefore different I-front expansion histories. We recall 
that the true solution of the problem is unknown. Given 
the non-linearity of the system, the differences in detail of 
the initial conditions and cooling rates, and the fact that we 
compare fundamentally different RT schemes, the agreement 
that we obtain is actually very good. 

To illustrate the spatial distribution of the ionised frac- 
tion and the temperature of the gas we show in Figure 
fT5l slices through the simulation volume at z = O.T^^box 
(through the largest group) at three different times t = 
0.05, 0.2 and 0.4 Myr. The upper row shows contours of the 
neutral fraction plotted over a density field, the second row 
shows a map of the neutral fraction and the third row shows 
a map of the temperature of the gas. The dense regions trap 
the I-front and thus produce sharp gradients of the radiation 
density. In the under-dense regions ionisation is more effec- 
tive and the I-front is extended. Even though the ionised re- 
gions are mostly uniform, traces of the dense structures that 
are less ionised can be seen near the front-trapping points. 
As shown in the contour maps, the I-fronts are broader in 
the low density regions and thiner in the high density re- 
gions. The temperature in the ionised regions reaches several 
10^ K and is uniform. It remains unchanged outside, where 
no photons are present. 

We further compare the volume faction of the temper- 
ature and the ionised fraction from GADGET, CRASH and 
FTTE in Figures [T6l and [T71 at three different times t = 0.05, 
0.2, and 0.4 Myr. The temperature volume fractions do not 
match particularly well due to the different photoheating 
mechanisms that GADGET, CRASH and FTTE use, but they 
find a similar maximum temperature. The volume fractions 
of the ionised fraction for GADGET, CRASH and FTTE have 
similar shapes. GADGET, in contrast to the other codes, pro- 
duces less intermediately ionised gas. However, overall the 
histograms are in a reasonably good agreement with each 



other, suggesting that our moment-based scheme is quite ca- 
pable in describing the reionisation process and produces re- 
sults of similar accuracy as other established radiative trans- 
fer codes. 



6 SUMMARY AND CONCLUSIONS 

We have presented a novel method for solving the radiative 
transfer equations within SPH, which is based on moments 
of the radiative transfer equation that are closed with a vari- 
able Eddington tensor. The radiation transport effectively 
becomes an anisotropic diffusion problem in this formula- 
tion. We have developed a new discretization scheme for 
anisotropic diffusion in SPH together with an implicit time 
integration method which for the first time allows a calcu- 
lation of such anisotropic diffusion in SPH. Together with 
a scheme to estimate Eddington tensors based on the opti- 
cally thin approximation, this yields a very fast approximate 
treatment of radiative transfer that can be used in dynami- 
cal SPH calculations. 

We have implemented our method into the cosmolog- 
ical simulation code GADGET-3 and presented several test 
problems where we varied the initial conditions and differ- 
ent numerical parameters to investigate the accuracy of the 
method. Our test results agree in general very well with an- 
alytical predictions and data from other simulations, except 
that the long-term evolution of sharp geometric shadows 
is clearly not followed accurately. While this clearly limits 
the range of applicability of the method, we expect that 
the method can still provide reasonably accurate results for 
problems where shadowing is comparatively unimportant, 
such as cosmological reionisation, where the SPH-based vari- 
able Eddington tensor approach can be competitive with 
other techniques. However, our method has two important 
strengths not shared by most other techniques: It can easily 
cope with an arbitrary number of sources since its speed is 
essentially independent of the number of sources, and fur- 
thermore, it is fast enough to be included into a cosmolog- 
ical simulation code where radiative transfer is calculated 
on-the-fly together with the ordinary dynamics. This is es- 
pecially promising for future calculations of galaxy forma- 
tion and reionisation that we want to carry out with our 
new code. 

The epoch of reionisation is an important and extremely 
interesting process in the evolution of the Universe. There 
is already very important observational evidence that con- 
strains reionisation, but few precise statements about the 
onset, duration, and end of reionisation can be made at this 
point. The situation will likely change in the near future 
through upcoming observations, ranging from new CMB ob- 
servations with PLANCK, to 21cm tomography at high red- 
shift with radio telescopes such as LOFAR. In the meantime, 
we heavily rely on cosmological simulations to advance our 
understanding of the reionisation process. 

However, the simulations run so far leave many ques- 
tions still unanswered, and their lack of self-consistency with 
the actual dynamics and the rapidly evolving source popu- 
lation may have introduced sizable inaccuracies. It is there- 
fore highly desirable to have new numerical techniques that 
do not rely on solving the radiative transfer equation in a 
post-processing mode by treating only individual snapshot 
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Figure 15. The upper row shows contours of neutral fractions equal to 0.01,0.5 and 0.9 in a slice through the simulation volume at 
z = O.T^box: through the largest group. The snapshots are taken at times t = 0.05, 0.2 and 0.4 Myr (left to right). The background shows 
a slice of the density distribution. The ionised regions expand with time as the I-front is trapped at high density regions and extends into 
low density regions. The second row shows the neutral fraction in the same slice. The ionised regions are uniform with some substructure 
visible near the front-trapping regions, where the I-front is not as diffuse as in low density regions. The third row shows the temperature 
distribution in the slice. The temperature in the ionised regions reaches several lO^K and remains uniform outside these regions. 






Figure 16. Volume fraction of the ionised fraction at three differ ent times t = 0. 05, 0.2, and 0.4 Myr (left to right). Results from 
GADGET are compared with results from CRASH and FTTE from llHev et"all (l2006l V All codes match in the shape of the histograms, 
but GADGET gives a lower intermediately ionised fraction. 



files from simulations. Rather, accurate calculations should 
account for radiative transfer 'on-the-fly', so that the gas 
properties are affected simultaneously by gravity and radi- 
ation. Our new method promises to be an important step 
in this direction. In future work, one of our most important 
goals is therefore to carry out high-resolution simulations of 



structure formation where the build-up of stellar and quasar 
sources in galaxies is followed self-consistently with the ra- 
diation field, allowing us to make more accurate predictions 
for the temporal evolution of reionisation, and the topology 
of reionised regions as a function of time. 

All tests problems we have presented in this thesis agree 
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Figure 17. Volume fraction of the temperature at thre e different t i mes t = 0.05, 0.2 and 0.4 Myr (left to right). Results from GADGET are 
compared with results from CRASH and FTTE from llliev et al.l (|2006l V All codes produce different histograms. FTTE and GADGET 
agree better with each other at later evolution times. CRASH produces higher temperatures due to its use of a different spectral 
distribution. 



well with theoretical predictions or results obtained with 
other radiative transfer codes. We should be able to obtain a 
realistic and accurate temperature evolution of the Universe 
during reionisation, which is important for setting the 'cos- 
mic equation of state' that regulates the absorption seen as 
Lyman-a forest in the spectra of distant quasars. Finally, we 
plan to include the photoionisation of other elements besides 
hydrogen, most importantly of helium. Helium reionisation 
probably happened sometime at redshift z ^ 2 — 4, where 
it may have left a sizable imprint in the temperature evolu- 
tion of the intergalactic medium. Surprisingly, recent obser- 
vations suggest that th e temperature-dens ity relation of the 
IGM may be inverted (|Bolton et al.ll2008l ). which could be 
caused by radiative transfer effects related to helium reioni- 
sation. Whether this is indeed possible can only be clarified 
with simulations. Studying this question with our new meth- 
ods would therefore be particularly timely. 
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APPENDIX A: COOLING RATES 

Here we present all cooling rates that we have used in our 
simulations with temperature evolution ( Section 15. 2p . The 
rates have been obtained from ICenI and are given in 

erg cm~^ s~"^. 



(i) Recombination cooling rate 



: 8. 7x10-2'^ I 



-V 



(los; 

(ii) Collisional ionisation cooling rate 

/ -157809.1 
105 ) V T 



neum , (Al) 



1. 27x10"^^ Vt 1 - 



neuiii, (A2) 



(iii) Collisional excitation cooling rate 

-1 



: 7.5x10-^^ 1 



118348 \ 




nefim, (A3) 



(iv) Bremsstrahlung cooling rate 
Ab = 1.42 X IQ-'^^gsVTnenm, 

where ^ff = 1.3 is the Gaunt factor. 



(A4) 



APPENDIX B: 
DESCENT 



METHOD OF STEEPEST 



In this Appendix, we give a derivation of the CG technique 
for solving linear problems, which we employ in our implicit 
time integration scheme of the anisotropic diffusion prob- 
lem. As the CG scheme is closely related to the method of 
steepest decent, we start with an explanation of this more 
general technique, and then specialize to the CG method. 

The method of steepest descent is a scheme to solve a 
linear system of equations given by Ax = b. The idea is to 
obtain the solution as the minimum of the quadratic form 



/(x) = ^x^^x - b^x, 

such that 
/'(x) 

This equation reduces to 
/'(x) = ^x-b 



-A x+-^x- 



(Bl) 



(B2) 



(B3) 



if A is symmetric and positive definite, i.e. if x^Ax > for 
all X / 0. 

We now consider an iteration scheme that tries to find 
the solution x. As we take steps, we choose the direction of 
the next step in the direction in which the quadratic form / 
decreases most rapidly, which is in the direction opposite to 
the gradient /^(x). Therefore, the next step should be pro- 
portional to — /^(x(^)) = b — Ax(^). Here the index i denotes 
the number of the step we take towards the correct value 
of X. Let us denote the difference between the numerical 



and the exact solution as e 



X, and the residual 



as r(i) = b — Ax(^) — — /^(x(i)). Therefore, the next step 
taken is given by X(i) = X(o) + Oi^(o) ■ The optimum value 
of a is chosen such that the directional derivative ^/(x(i)) 
equals 0, i.e. the vectors /^(x(i)) and r(o) should be chosen 
orthogonal: 



— /(x(i)) = /(x(i))^— X(i) = /(x(i))^r(o) = 



da 



(B4) 
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We further notice that /^(x(i)) = — and therefore 

rfi)r(o) = 

(b- Ax(i))^r(o) = 

(b- A(x(o) +ar(o)))'^r(o) = 

(b - Ax(o))'^r(o) - a(Ar(o))'^r(o) = 

(b - Ax(o))^r(o) = a(Ar(o))^r(o) 

rfo)r(o) = arfo)(Ar(o)) 

^fo)^(o) 
rfo)Ar(o) 

Finally, putting it all together, the method of steepest de- 
scent is as follows 

r(i) = b-Ax(i) (B5) 

T 

X(i+i) = X(i) +a(i)r(i) (B7) 
r(i+i) = r(i) - Q;(i)Ar(i) (B8) 

APPENDIX C: CONJUGATE GRADIENT 
METHOD 

In the method of steepest descent the value of x is deter- 
mined via successively adding the search directions r(^) . Let 
us define the set of directions {d(4)} as the search directions 
for the CG method, such that X(^+i) = X(^) + a(^)d(^) and 
r(i+i) = r(^) — a(^i)Ad(^iy We further require that the vectors 
{d(i)} are A-conjugate, i.e. d(i)Ad(j) = 0, which means 

i 

d(i+i) = + ^ ftfcd(jt). (CI) 

A;=0 

Using Gram- Schmidt orthogonalization, the coefficients 



Pik are 


found to be 






Pik = < 


1 •^(i)^(^) 


for z = /c + 1 
for z > /c + 1 


(C2) 


Thus 










r(i_i)r(i-i) 




(C3) 


Therefore, the CG method can 


be summarized 


as follows: 


r(o) 


= b - Ax(o) 




(C4) 


d(o) 


= b - Ax(o) 




(C5) 




df.)Ad(i) 




(C6) 




= X(i) + Q;(i)d(i) 




(C7) 


r(i+i) 


= r(i) - Q;(i)Ad(i) 




(C8) 








(C9) 


d(i+i) 


= ^(i+i) ~ Ai+i)Ad( 




(CIO) 



The interesting feature of the CG method is that each 
subsequent correction to the solution vector is orthogonal 



to all previous ones, while at the same time it points into 
the direction where the error in the solution decreases most 
quickly. This normally produces a comparatively rapid con- 
vergence of the scheme. 
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